The present invention relates to processing of seismic data representative of subsurface features in the earth and, more particularly, to improved methods of processing seismic data using high resolution Radon transformations to remove unwanted noise from meaningful reflection signals.
Seismic surveys are one of the most important techniques for discovering the presence of oil and gas deposits. If the data are properly processed and interpreted, a seismic survey can give geologists a picture of subsurface geological features, so that they may better identify those features capable of holding oil and gas. Drilling is extremely expensive, and ever more so as easily tapped reservoirs are exhausted and new reservoirs are harder to reach. Having an accurate picture of an area's subsurface features can increase the odds of hitting an economically recoverable reserve and decrease the odds of wasting money and effort on a nonproductive well.
The principle behind seismology is deceptively simple. As seismic waves travel through the earth, portions of that energy are reflected back to the surface as the energy waves traverse different geological layers. Those seismic echoes or reflections give valuable information about the depth and arrangement of the formations, some of which hopefully contain oil or gas deposits.
A seismic survey is conducted by deploying an array of energy sources and an array of sensors or receivers in an area of interest. Typically, dynamite charges are used as sources for land surveys, and air guns are used for marine surveys. The sources are discharged in a predetermined sequence, sending seismic energy waves into the earth. The reflections from those energy waves or “signals” then are detected by the array of sensors. Each sensor records the amplitude of incoming signals over time at that particular location. Since the physical location of the sources and receivers is known, the time it takes for a reflection signal to travel from a source to a sensor is directly related to the depth of the formation that caused the reflection. Thus, the amplitude data from the array of sensors can be analyzed to determine the size and location of potential deposits.
This analysis typically starts by organizing the data from the array of sensors into common geometry gathers. That is, data from a number of sensors that share a common geometry are analyzed together. A gather will provide information about a particular spot or profile in the area being surveyed. Ultimately, the data will be organized into many different gathers and processed before the analysis is completed and the entire survey area mapped.
The types of gathers typically used include: common midpoint, where the sensors and their respective sources share a common midpoint; common source, where all the sensors share a common source; common offset, where all the sensors and their respective sources have the same separation or “offset”; and common receiver, where a number of sources share a common receiver. Common midpoint gathers are the most common gather today because they allow the measurement of a single point on a reflective subsurface feature from multiple source-receiver pairs, thus increasing the accuracy of the depth calculated for that feature.
The data in a gather are typically recorded or first assembled in the offset-time domain. That is, the amplitude data recorded at each of the receivers in the gather are assembled or displayed together as a function of offset, i.e., the distance of the receiver from a reference point, and as a function of time. The time required for a given signal to reach and be detected by successive receivers is a function of its velocity and the distance traveled. Those functions are referred to as kinematic travel time trajectories. Thus, at least in theory, when the gathered data are displayed in the offset-time domain, or “X-T” domain, the amplitude peaks corresponding to reflection signals detected at the gathered sensors should align into patterns that mirror the kinematic travel time trajectories. It is from those trajectories that one ultimately may determine an estimate of the depths at which formations exist.
A number of factors, however, make the practice of seismology and, especially, the interpretation of seismic data much more complicated than its basic principles. First, the reflected signals that indicate the presence of geological strata typically are mixed with a variety of noise.
The most meaningful signals are the so-called primary reflection signals, those signals that travel down to the reflective surface and then back up to a receiver. When a source is discharged, however, a portion of the signal travels directly to receivers without reflecting off of any subsurface features. In addition, a signal may bounce off of a subsurface feature, bounce off the surface, and then bounce off the same or another subsurface feature, one or more times, creating so-called multiple reflection signals. Other portions of the signal turn into noise as part of ground roll, refractions, and unresolvable scattered events. Some noise, both random and coherent, is generated by natural and man-made events outside the control of the survey.
All of this noise is occurring simultaneously with the reflection signals that indicate subsurface features. Thus, the noise and reflection signals tend to overlap when the survey data are displayed in X-T space. The overlap can mask primary reflection signals and make it difficult or impossible to identify patterns in the display upon which inferences about subsurface geology may be drawn. Accordingly, various mathematical methods have been developed to process seismic data in such a way that noise is separated from primary reflection signals.
Many such methods seek to achieve a separation of signal and noise by transforming the data from the X-T domain to other domains. In other domains, such as the frequency-wavenumber (F-K) domain or the time-slowness (tau-P), there is less overlap between the signal and noise data. Once the data are transformed, various mathematical filters are applied to the transformed data to eliminate as much of the noise as possible and, thereby, to enhance the primary reflection signals. The data then are inverse transformed back into the offset-time domain for interpretation or further processing.
For example, so-called Radon filters are commonly used to attenuate or remove multiple reflection signals. Such methods rely on Radon transformation equations to transform data from the offset-time (X-T) domain to the time-slowness (tau-P) where it can be filtered. More specifically, the X-T data are transformed along kinematic travel time trajectories having constant velocities and slownesses, where slowness p is defined as reciprocal velocity (or p=1/v).
Such prior art Radon methods, however, typically first process the data to compensate for the increase in travel time as sensors are further removed from the source. This step is referred to as normal moveout or “NMO” correction. It is designed to eliminate the differences in time that exist between the primary reflection signals recorded at close-in receivers, i.e., at near offsets, and those recorded at remote receivers, i.e., at far offsets. Primary signals, after NMO correction, generally will transform into the tau-P domain at or near zero slowness. Thus, a mute filter may be defined and applied in the tau-P domain. The filter mutes, i.e., excludes all data, including the transformed primary signals, below a defined slowness value pmute.
The data that remain after applying the mute filter contains a substantial portion of the signals corresponding to multiple reflections. That unmuted data are then transformed back into offset-time space and are subtracted from the original data in the gather. The subtraction process removes the multiple reflection signals from the data gather, leaving the primary reflection signals more readily apparent and easier to interpret.
It will be appreciated, however, that in such prior art Radon filters, noise and multiple reflection signals recorded by receivers close to the gather reference point (“near trace multiples”) are not as effectively separated from the primary reflections. The lack of separation in large part is an artifact of the NMO correction performed prior to transformation. Because NMO correction tends to eliminate offset time differences in primary signals, primary signals and near trace multiple signals both transform at or near zero slowness in the tau-P domain. When the mute filter is applied, therefore, the near trace multiples are muted along with the primary signal. Since they are muted, near trace multiples are not subtracted from and remain in the original data gather as noise that can mask primary reflection signals.
Radon filters have been developed for use in connection with common source, common receiver, common offset, and common midpoint gathers. They include those based on linear slant-stack, parabolic, and hyperbolic kinematic travel time trajectories. The general case forward transformation equation used in Radon filtration processes, R(p,τ)[d(x,t)], is set forth below:
            u      ⁡              (                  p          ,          τ                )              =                  ∫                  -          ∞                ∞            ⁢                        ⅆ          x                ⁢                              ∫                          -              ∞                        ∞                    ⁢                                    ⅆ              t                        ⁢                                                  ⁢                          d              ⁡                              (                                  x                  ,                  t                                )                                      ⁢                          δ              ⁡                              [                                  f                  ⁡                                      (                                          t                      ,                      x                      ,                      τ                      ,                      p                                        )                                                  ]                                                              (          forward      ⁢                          ⁢      transformation        )  where                u(p,τ)=transform coefficient at slowness p and zero-offset time τ        d(x,t)=measured seismogram at offset x and two-way time t        p=slowness        t=two-way travel time        τ=two-way intercept time at p=0        x=offset        δ=Dirac delta function        ƒ(t,x,τ,p)=forward transform function        
The forward transform function for linear slant stack kinematic travel time trajectories is as follows:
            f      ⁡              (                  t          ,          x          ,          τ          ,          p                )              =          t      -      τ      -              p        ⁢                                  ⁢        x              where                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                    ⁢                      δ            ⁡                          (                              t                -                τ                -                                  p                  ⁢                                                                          ⁢                  x                                            )                                                                                =                        ⁢            1                    ,                                    when              ⁢                                                          ⁢              t                        =                          τ              +                              p                ⁢                                                                  ⁢                x                                              ,          and                                                          =                        ⁢            0                    ,                      elsewhere            .                              Thus, the forward linear slant stack Radon transformation equation becomes
      u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢                          ⁢              d        ⁡                  (                      x            ,                          τ              +                              p                ⁢                                                                  ⁢                x                                              )                    
The forward transform function for parabolic kinematic trajectories is as follows:
            f      ⁡              (                  t          ,          x          ,          τ          ,          p                )              =          t      -      τ      -              p        ⁢                                  ⁢                  x          2                      where                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                    ⁢                      δ            ⁡                          (                              t                -                τ                -                                  p                  ⁢                                                                          ⁢                                      x                    2                                                              )                                                                                =                        ⁢            1                    ,                                    when              ⁢                                                          ⁢              t                        =                          τ              +                              p                ⁢                                                                  ⁢                                  x                  2                                                              ,          and                                                          =                        ⁢            0                    ,                      elsewhere            .                              Thus, the forward parabolic Radon transformation equation becomes
      u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢                          ⁢              d        ⁡                  (                      x            ,                          τ              +                              p                ⁢                                                                  ⁢                                  x                  2                                                              )                    
The forward transform function for hyperbolic kinematic travel time trajectories is as follows:
            f      ⁡              (                  t          ,          x          ,          τ          ,          p                )              =          t      -                                    τ            2                    -                                    p              2                        ⁢                          x              2                                            where                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                    ⁢                      δ            ⁡                          (                              t                -                                                                            τ                      2                                        +                                                                  p                        2                                            ⁢                                              x                        2                                                                                                        )                                                                                =                        ⁢            1                    ,                                    when              ⁢                                                          ⁢              t                        =                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                ,          and                                                          =                        ⁢            0                    ,                      elsewhere            .                              Thus, the forward hyperbolic Radon transformation equation becomes
      u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢                          ⁢              d        ⁡                  (                      x            ,                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
A general case inverse Radon transformation equation is set forth below:
            d      ⁡              (                  x          ,          t                )              =                  ∫                  -          ∞                ∞            ⁢                        ⅆ          p                ⁢                              ∫                          -              ∞                        ∞                    ⁢                                    ⅆ                              τρ                ⁡                                  (                  τ                  )                                                      *                          u              ⁡                              (                                  p                  ,                  τ                                )                                      ⁢                          δ              ⁡                              [                                  g                  ⁡                                      (                                          t                      ,                      x                      ,                      τ                      ,                      p                                        )                                                  ]                                                              (          inverse      ⁢                          ⁢      transformation        )  where                g(t,x,τ,p)=inverse transform function, and        ρ(τ)*=convolution of rho filter.        
The inverse transform function for linear slant stack kinematic trajectories is as follows:g(t,x,τ,p)=τ−t+pxThus, when τ=t−px, the inverse linear slant stack Radon transformation equation becomes
      d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢                          ⁢              ρ        ⁡                  (          τ          )                    *              u        ⁡                  (                      p            ,                          t              -                              p                ⁢                                                                  ⁢                x                                              )                    
The inverse transform function for parabolic trajectories is as follows:g(t,x,τ,p)=τ−t+px2Thus, when τ=t−px2, the inverse parabolic Radon transformation equation becomes
      d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢                          ⁢              ρ        ⁡                  (          τ          )                    *              u        ⁡                  (                      p            ,                          t              -                              p                ⁢                                                                  ⁢                                  x                  2                                                              )                    
The inverse transform function for hyperbolic trajectories is as follows:g(t,x,τ,p)=τ−√{square root over (t2−p2x2)}Thus, when τ=√{square root over (t2−p2x2)}, the inverse hyperbolic Radon transformation equation becomes
      d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢                          ⁢              ρ        ⁡                  (          τ          )                    *              u        ⁡                  (                      p            ,                                                            t                  2                                -                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
The choice of which form of Radon transformation preferably is guided by the travel time trajectory at which signals of interest in the data are recorded. Common midpoint gathers, because they offer greater accuracy by measuring a single point from multiple source-receiver pairs, are preferred over other types of gathers. Primary reflection signals in a common midpoint gather generally will have hyperbolic kinematic trajectories. Thus, it would be preferable to use hyperbolic Radon transforms.
To date, however, Radon transformations based on linear slant stack and parabolic trajectories have been more commonly used. The transform function for hyperbolic trajectories contains a square root whereas those for linear slant stack and parabolic transform functions do not. Thus, the computational intensity of hyperbolic Radon transformations has in large part discouraged their use in seismic processing.
It has not been practical to accommodate the added complexity of hyperbolic Radon transformations because the computational requirements of conventional processes are already substantial. NMO correction involves a least-mean-squares (“LMS”) energy minimization calculation. Forward and inverse Radon transforms also are not exact inverses of each other. Accordingly, an additional LMS calculation is often used in the transformation. Those calculations in general and, especially LMS analyses, are complex and require substantial computing time, even on advanced computer systems.
Moreover, a typical seismic survey may involve hundreds of sources and receivers, thus generating tremendous quantities of raw data. The data may be subjected to thousands of different data gathers. Each gather is subjected not only to filtering processes as described above, but in all likelihood to many other enhancement processes as well. For example, data are typically processed to compensate for the diminution in amplitude as a signal travels through the earth (“amplitude balancing”). Then, once the individual gathers have been filtered, they must be “stacked”, or compiled with other gathers, and subjected to further processing in order to make inferences about subsurface formations. Seismic processing by prior art Radon methods, therefore, requires substantial and costly computational resources, and there is a continuing need for more efficient and economical processing methods.
It is understood, of course, that the transformation, NMO correction, and various other mathematical processes used in seismic processing do not necessarily operate on all possible data points in the gather. Instead, and more typically, those processes may simply sample the data. For example, the transform functions in Radon transformations are based on variables t, x, τ, and p. The transformations, however, will not necessarily be performed at all possible values for each variable. Instead, the data will be sampled a specified number of times, where the number of samples for a particular variable may be referred to as the index for the variable set. For example, k, l, m, and j may be defined as the number of samples in, or the index of, respectively, the time (t), offset (x), tau (τ), and slowness (p) sets. The samples will be taken at specified intervals, Δt, Δx, Δτ, and Δp. Those intervals are referred to as sampling variables.
The magnitude of the sampling variables, the variable domain and all other factors being equal, will determine the number of operations or samples that will be performed in the transformation and the accuracy of the transformation. As the sampling variables become finer, i.e., smaller, the number of samples increases and so does the accuracy. By making the sampling variables larger, or coarser, the number of samples is reduced, but the accuracy of the transformation is reduced as well.
In large part because of the computational intensity of those processes, however, the sampling variables used in prior art Radon processes can be relatively coarse. In NMO correction, for example, industry typically samples at Δt values several time greater than that at which the field data were recorded, i.e., in the range of 20-40 milliseconds, and samples at v values of from about 15 to about 120 m/sec. To the extent that the sampling variables are greater than the data acquisition sample rates, data are effectively discarded and the accuracy of the process is diminished. In the Radon transformation itself, prior art methods typically set Δt, Δτ, and Δx at the corresponding sample rates at which the data were recorded, i.e., in the range of 2 to 4 milliseconds and 25 to 400 meters or more, respectively. Δp is relatively coarse, however, typically being set at values in the range of 4-24 μsec/m. Also, the index j of the slowness (p) set usually is equal to the fold of the offset data, i.e., the number of source-receiver pairs in the survey, which typically ranges from about 20 to about 120.
Accordingly, prior art methods may rely on additional computations in an attempt to compensate for the loss of accuracy and resolution inherent in coarse sampling variables. For example, the LMS analysis applied in NMO correction and prior art Radon transformations are believed to obviate the need for finer sampling variables. Certain prior art processes utilize trace interpolation, thereby nominally reducing x values, but the additional “data” are approximations of what would be recorded at particular offsets between the offsets where receivers were actually located. Relative inaccuracy in processing individual data gathers also is believed to be acceptable because the data from a large number of gathers ultimately are combined for analysis, or “stacked”. Nevertheless, there also is a continuing need for more accurate processing of seismic data and processes having finer resolutions.
The velocity at which primary reflection signals are traveling, what is referred to as the stacking velocity, is used in various analytical methods that are applied to seismic data. For example, it is used in determining the depth and lithology or sediment type of geological formations in the survey area. It also is used various seismic attribute analyses. A stacking velocity function may be determined by what is referred to as a semblance analysis. Which such methods have provided reasonable estimations of the stacking velocity function for a data set, given the many applications based thereon, it is desirable to define the stacking velocity function as precisely as possible.
An object of this invention, therefore, is to provide a method for processing seismic data that more effectively removes unwanted noise from meaningful reflection signals.
It also is an object to provide such methods based on hyperbolic Radon transformations.
Another object of this invention is to provide methods for removal of noise from seismic data that are comparatively simple and require relatively less computing time and resources.
Yet another object is to provide methods for more accurately defining the stacking velocity function for a data set.
It is a further object of this invention to provide such methods wherein all of the above-mentioned advantages are realized.
Those and other objects and advantages of the invention will be apparent to those skilled in the art upon reading the following detailed description and upon reference to the drawings.